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NUMERICAL SOLUTION OF THE EQUATIONS 
FOR COMPRESSIBLE LAMINAR, TRANSITIONAL, AND 
TURBULENT BOUNDARY LAYERS AND COMPARISONS 
WITH EXPERIMENTAL DATA* 

By Julius E. Harris 
Langley Research Center 

SUMMARY 

A numerical method for solving the system of equations which govern the mean flow 
properties of laminar, transitional, and turbulent compressible boundary layers for either 
planar or axisymmetric flows is presented. The turbulent boundary layer is treated by 
a two-layer concept with appropriate eddy viscosity models used for each layer to replace 
the Reynolds stress term. A specifiable static turbulent Prandtl number relates the 
turbulent heat flux term to the Reynolds stress. The mean properties in the transitional 
boundary layer are calculated by multiplying the eddy viscosity by an intermittency func- 
tion based on the statistical production and growth of the turbulent spots. The numerical 
method used to solve the system of equations is a three-point implicit finite -difference 
scheme. The momentum and energy equations are solved simultaneously without iteration. 

A number of test cases are compared with experimental data for supersonic and 
hypersonic flows over planar and axisymmetric geometries. These test cases include 
laminar, transitional, and turbulent boundary-layer flows with both favorable and unfavor- 
able pressure -gradient histories. Mass flux at the wall and transverse curvature effects 
are considered. The results show that the system of equations and the numerical tech- 
nique provide accurate predictions for laminar, transitional, and turbulent compressible 
boundary-layer flows. 


*The information presented herein is based in part upon a thesis entitled "Numeri- 
cal Solution of the Compressible Laminar, Transitional, and Turbulent Boundary Layer 
Equations With Comparisons to Experimental Data," submitted in partial fulfillment of 
the requirements for the degree of Doctor of Philosophy in Aerospace Engineering, 
Virginia Polytechnic Institute and State University, Blacksburg, Virginia, May 1970. 



INTRODUCTION 


The boundary-layer concept first introduced by Prandtl (ref. 1) in 1904 divides the 
flow field over an arbitrary surface into two distinct regions; an inviscid outer region in 
which solutions to the Euler equations describe the flow-field characteristics, and a 
viscous inner region where the classical boundary-layer equations apply. The boundary- 
layer region may be further categorized as to type, namely, laminar, transitional, and 
turbulent. 

The laminar boundary layer has received considerable attention over the past 
60 years. Reviews of early methods are presented in references 2 to 4. A review of 
similar and local similarity solutions is presented in reference 5. In the early part of 
the past decade, the complete nonsimilar laminar equations were solved to a high degree 
of accuracy by finite -difference techniques. (See Blottner, ref. 6.) 

The mean flow within the transition region has not been studied as extensively as 
either the location of transition or the characteristics of the fully developed turbulent 
boundary layer. There have been many experimental studies where the heat transfer at 
the wall was measured, but these studies reveal little of the flow structure away from 
the wall. A few experimental studies have been made in which the instantaneous behavior 
of the transitional boundary layer was reported for both incompressible (refs. 7 to 9) and 
compressible flows (ref. 10) but more detailed work is still required. The available 
transitional region data, although limited in extent, does allow workable models of the 
mean flow structure in the transition region to be formulated and applied tentatively to 
compressible flow systems. These crude models are generally based on the transport 
models for fully developed turbulent flows modified by some intermittency distribution. 
Reviews of transition and the transitional region have recently been presented by 
Savulescu (ref. 11) and Morkovin (ref. 12). (See, also, ref. 13.) 

Compressible turbulent boundary-layer flows have received accelerated study over 
the past decade. At first most of the work was experimental, the main objective being 
development of empirical or semiempirical correlation techniques. Little effort was 
devoted to obtaining numerical solutions of the equations for turbulent boundary layers 
until a few years ago. The principal difficulties were associated with modeling the 
turbulent transport terms as well as developing techniques for obtaining solutions on 
existing digital computer systems. Even today, because of the limited understanding of 
these turbulent transport processes, completely general solutions of the mean turbulent 
boundary -layer equations are not possible. However, by modeling the turbulent transport 
terms through eddy viscosity or mixing-length concepts, it is possible to solve the sys- 
tem of equations directly. Reviews of recent analytical advances are contained in 
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reference 14 for incompressible flow. (See ref. 15 for compilation of incompressible 
data.) Similar material for compressible flows is presented in references 16 and 17. 

A number of papers have been presented over the past decade in which attempts 
have been made to connect the three boundary-layer flow regions (laminar, transitional, 
and turbulent) so that one system of equations could be used to describe compressible 
boundary-layer flows. Persh (ref. 18) appears to be one of the first to develop reason- 
ably accurate solutions for the transitional region. Persh used an empirical correlation 
of velocity profile data together with the streamwise momentum equation to predict the 
characteristics of transitional flow; although limited in application, it did present a 
method with which the laminar and turbulent regions could be tied together. Weber 
(ref. 19) used mixing-length theory and an intermittency factor based on the probability 
equation of Emmons (ref. 20) to calculate velocity profiles which were in fair agreement 
with existing experimental data. (See also Masaki and Yakura (ref. 21).) Martellucci, 

Rie, and Sontowski (ref. 22) were among the first to define an effective viscosity model for 
laminar, transitional, and turbulent compressible boundary-layer flows. (See also, 

Harris (ref. 23), Adams (ref. 24), and Kuhn (ref. 25)). Fish and McDonald (ref. 26) have 
developed a solution technique in which the mixing length at each streamwise station is 
governed by an integral solution of the turbulent kinetic energy equation with the free- 
stream turbulence intensity used as a boundary condition. The most significant feature 
of the Fish-McDonald technique is that the transition location and the extent of the transi- 
tional flow are not required as inputs to the solution. The method has currently been 
applied only to incompressible boundary-layer flows. (See also, Glushko (ref. 27) and 
Donaldson (ref. 28)). 

The present paper presents a system of equations that are applicable to laminar, 
transitional, and turbulent compressible boundary layers and a numerical technique for 
obtaining accurate solutions of the system for either planar or axisymmetric perfect gas 
flows. The numerical technique is similar to that first developed by Fliigge-Lotz and 
Blottner (ref. 29) and later improved by Davis and Fliigge-Lotz (ref. 30). The momentum 
and energy equations remain coupled and are solved simultaneously in the transformed 
plane without iteration. (Difference relations and coefficients for the difference equa- 
tions are developed in appendixes A and B, respectively.) The procedure's very effi- 
cient for parabolic systems where only two governing equations are required; however, 
as the number of equations increases beyond two, the method becomes increasingly 
inefficient because of the digital computer storage requirements. (See ref. 6.) The 
numerical technique yields accurate results for compressible laminar, transitional, and 
turbulent boundary layers with pressure gradients, heat transfer, and mass transfer at 
the wall. The method treats the fully developed turbulent region by replacing the Reynolds 
stress terms with an eddy viscosity model. A specifiable static turbulent Prandtl num- 
ber function relates the turbulent flux of heat to the eddy viscosity model. The mean 
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properties of the transitional boundary layer are calculated by multiplying the eddy vis- 
cosity by an intermittency function based on the statistical production and growth of the 
turbulent spots. The eddy viscosity model is based upon existing experimental data. 


SYMBOLS 


A damping function (eq. (38)) 

A (1) b* 1 ) c (1 > 

A n ,t5 n ’Si ’ 

D (l) E (D f (D G d) 

u n j^n > f n >u n 


coefficients in difference equation (66) and defined by equa- 
tions (B3) to (B9) 


a ( 2 ) b ( 2 ) c ( 2 ) 

A n ’Si ’ c n ’ 

D (2) e (2) f < 2) g < 2 > 
n ’n ’ n ’n 


coefficients in difference equation (67) and defined by equa- 
tions (BIO) to (B16) 


Cf 


skin-friction coefficient, 


'w 


(a pu )e 


c ml’ C ml 

C P 


c 


®ml)^ml 


Ey 

F 

/V 

F 


Fml 

F m2 

Fy 


defined in equations (B45) and (B46) 

specific heat at constant pressure 

length defined in figure 13(a) 

defined in equations (B36) and (B37) 

defined in equation (B39) 

velocity ratio, — 
u e 

(pv) w 

mass injection parameter, 

(pu) e 

defined in equation (B29) 
defined in equation (B32) 
defined in equation (B40) 


4 



G 


a typical quantity in boundary layer 


H 

h 15 h 2 ,h 3 , 


h 

i 

3 

K 


Kl 


k 2 

K l 

k t 

L 

l 

l 

M 


m 

N 

N Pr 

N Pr,t 


a typical quantity in boundary layer (see appendix A) 

. coefficients defined by equations (B17) to (B28) 

heat -transfer coefficient 
index used in grid-point notation (eq. (64)) 
flow index; j = 0 planar flow, j = 1 axisymmetric flow 
grid-point -spacing parameter (eq. (64)) 
constant in eddy viscosity model (eq. (36)) 
constant in eddy viscosity model (eq. (41)) 
thermal conductivity 
eddy conductivity (eq. (9)) 
reference length 

defined in equations (B34) and (B35) 
defined in equations (28) 
mixing length (eq. (36)) 

Mach number 
grid-point index (fig. 7) 

number of grid points at each x-station (fig. 7) 

C d M 

Prandtl number, — t — 

static turbulent Prandtl number (eq. (10)) 


N ! 

n 

P 

P 

Q 

q 

R 

R. 


St 


R, 


e,x 


R 


e > x t,i 


R, 


e,Ax t 


R .* 

e,6t 


R 


e,6 


R 


°°,d 


S 

T 


Stanton number, 


c p pu 


index defined in figure 7 


coefficient in correlation equation (55) 


pressure 


coefficient in correlation equation (55) 


heat -transfer rate 


coefficient in correlation equation (55) 


unit Reynolds number, — 

v e 

Reynolds number based on x, 


u e,x 

v e 


Reynolds number at transition, — x, . 

v e t ’ 1 

Up 

Reynolds number based on transition extent, —(x. t - x, ,\ 

’ v e \ t,f t,ij 

transition Reynolds number based on displacement thickness, 


Reynolds number based on momentum thickness, 1-^)0 




u d 

free-stream Reynolds number based on diameter, —2- 

^ on 


radial coordinate (fig. 1) 


recovery factor (eq. (75)) 


nose radius 


body radius (fig. 1) 


Sutherland viscosity constant (198.6° R (110.3° K)) 


static temperature 



defined in equations (B30) and (B33) 
defined in equation (B41) 




u 


u + 


U, 


transverse curvature term (eqs. (21)) 
velocity component in x-direction (fig. 1) 
law-of-wall coordinate, — 

U T 

friction velocity, \j r w/P 
V transformed normal velocity component (eq. (24)) 

V ml defined in equation (B31) 

v velocity component in y-direction 

v velocity component defined by equation (7) 

X 1 ,X 2 ,X 3 ,X 4 ,X 5 functions of grid-point spacing defined by equations (A4) to (A8) 

x boundary-layer coordinate tangent to surface 

Xj. ^ end of transition 

x, . beginning of transition 

functions of grid-point spacing defined by equations (A12) to (A17) 
Yr 7 ,Yg,Yg,YiQ functions defined by equations (70) 


y 

y + 

y m 


boundary-layer coordinate normal to surface 

yu 

law of wall coordinate, — - 


defined in figure 4 


body coordinate (fig. 1) 


a 


defined in equations (28) 



<3 

r 

y 

y 

Ax 

Ax t 

UiAv 

6 

6 * 

** 

®inc 

6 

6 

e av 

6 

7 
© 
e 
x 

v 
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defined in equations (28) 

streamwise intermittency distribution (eq. (57)) 
ratio of specific heats 

transverse intermittency distribution (eq. (44)) 
grid-point spacing in physical plane 
transition extent, f _ x t i 

grid-point spacing in transformed plane (see fig. 7) 
boundary-layer thickness 
displacement thickness 
incompressible displacement thickness 
eddy viscosity 

eddy viscosity function defined in equation (15) 
defined in equation (74) 

eddy viscosity function defined in equation (16) 

transformed normal boundary -layer coordinate 

static temperature ratio (eqs. (23)) 

momentum thickness 

defined in equation (59) 

molecular viscosity 

kinematic viscosity, 

r 



V 

l 


n > n l> n 2 

P 

T 

0 

X 


^max 

(Xmax) cr 

S2 


average kinematic viscosity 

transformed streamwise boundary-layer coordinate 

defined in equation (58) 

functions defined in equations (48) to (50) 

density 

shear stress 

angle between local tangent to surface and center line of body (see fig. 1) 

vorticity Reynolds number (eq. (54)) 

maximum local value of y (fig* 5) 

value of x max at transition (stability index) 

functional relation (see eq. (62)) 


Subscripts: 

e edge value or based on edge conditions 

i inner region of turbulent layer 

m mesh point in ^-direction (see fig. 7) 

max maximum value 

min minimum value 

n mesh point in ^-direction (see fig. 7) 

o outer region of turbulent layer 

s l sublayer edge 
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sp 


stagnation point 


t total condition 

u uniform or average value 

w wall value 

°° free stream 


Superscripts: 
j flow index 

* fluctuating component 

— time average value 


Other notation: 

TVC transverse curvature 

A coordinate used as a subscript means partial differential with respect to the 
coordinate. (See eq. (Al).) 

EQUATIONS FOR THE LAMINAR, TRANSITIONAL, AND 
TURBULENT COMPRESSIBLE BOUNDARY LAYER 

This section presents the governing equations for the compressible boundary layer 
together with the required boundary conditions. The eddy viscosity and eddy conductivity 
models used to represent the apparent turbulent shear and heat-flux terms appearing in 
the mean turbulent boundary-layer equations are discussed. 

Physical Plane 

Geometry and notation.- The orthogonal coordinate system is presented in figure 1. 
The boundary-layer coordinate system is denoted by X and Y which are tangent to and 
normal to the surface, respectively. The origin of both the boundary- layer coordinate 
axes system X,Y and the body coordinate axes system Z,R is located at the stagnation 
point for blunt bodies as shown in figure 1 or at the leading edge for sharp-tipped bodies. 
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Figure 1. - Coordinate system and notation. 


The velocity components u and v are oriented in the X- and Y-direction, respectively. 
Transverse curvature terms are retained because of their importance in the development 
of boundary -layer flows over slender bodies of revolution where the boundary -layer 
thickness may become of the order of the body radius r 0 . The angle 0 is the angle 
between the Z axis and local tangent evaluated at (x,0). The coordinates ^ and 
(xt £,0*) represent the location at which transition is initiated and completed, respectively.. 


Differential equations .- The flow of a compressible, viscous, heat conducting fluid 
is mathematically described by the continuity, Navier-Stokes, and energy equations 
together with an equation of state, a heat conductivity law, and a viscosity law. For flows 
at large Reynolds number, Prandtl (ref. 1) has shown that the Navier-Stokes and energy 
equations can be simplified to a form now recognized as the classical boundary-layer 
equations. These equations may be written as follows: 
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Continuity: 


i;( r3pu ) + -^( rjpv ) = 0 


( 1 ) 


Momentum: 


/ 9u „ 0u\ 

dp . 1 

_a_ 

j/ 

„ 9u \ 

U — + V — 

\ 9x 9yy 

n 

i 

ii 

9y 

r J 

1 9yj 


( 2 ) 


Energy: 


“ Jt( c p t ) + v #( c p t ) 


dp l 
= u — + — 
dx v i 


8 

8 

9y 

1 

•P 

•21 


9u 


(3) 


Osborne Reynolds (ref. 31) in 1883 was the first to observe and study the phenomena 
of transition from laminar to turbulent flow. Reynolds assumed that the instantaneous 
fluid velocity satisfied the Navier-Stokes equations and that the instantaneous velocity 
could be separated into mean and fluctuating components. The result of his early work 
was a set of Reynolds equations which differed from the Navier-Stokes equations only 
through the additional terms called the Reynolds stresses. These equations can be 
reduced through the boundary-layer approximations and written in terms of the mean 
variables as follows (ref. 32): 


Continuity: 





= 0 


Momentum: 


(4) 


Energy: 



(5) 


= „ dp 1 _8_T j JL (c T) 

dx + r j ay C p 9y Cp 



+ 



-,~7 8u 

" P U V ^ 


( 6 ) 
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The reduction of the Reynolds equations to the mean turbulent boundary-layer form 
presented in equations (4), (5), and (6) requires a number of stringent assumptions based 
on an order-of-magnitude analysis. These assumptions together with the various corre- 
lation terms which are neglected are discussed by Van Driest (ref. 32). However, for 
hypersonic flows where the boundary- layer thickness may increase rapidly and the density 
gradients are large, the order-of-magnitude analysis as present ed in re ference 32 may no 
longer be valid. For example, the Reynolds shear stress term (pv)'u' yields three cor- 
relation products when expanded; namely, p(u'v'), v(p'u'), and p'u’v'. The v(p'u') and 
p’u'v' correlations are generally neglected for supersonic flows since they are of lower 
order than the p(u’v') correlation; however, for hypersonic flows these two correlations 
may be of the same order or larger than the p(u'v') term and as such should be retained 
in the governing equations. (See Bushnell and Beckwith (ref. 33).) The term (pv)'T' 
appearing in the energy equation (ref. 32) may also be expanded into three correlation 
terms; namely, p(v'T'), v(p'T'), and p'v'T’. For supersonic flows the correlation 
terms v(p'T*) and p'v’T' are of lower order than p(v'T r ) and can be neglected; how- 
ever, for hypersonic flows these terms may again be of equal order and should then be 

r\ 

retained in the governing equations. Another correlation term ■^(p’v') which appears 
in the static energy equation ( ref. 32) is also neglected. For hypersonic flows this term 
may be of the same order as the Reynolds stress term retained in equation (5) and should 
then be retained in the system of equations. The actual effect of neglecting the various 
correlation terms can only be determined after well-documented, accurate experimental 
data become available for the hypersonic speed range. 

The mean turbulent equations are identical to those for laminar flows with the 
exception of the correlations of the turbulent fluctuating quantities which represent the 
apparent mass, shear, and heat-flux terms caused by the turbulence. The major problem 
encountered in calculating turbulent flows from this set of equations (eqs. (4) to (6)) is 
how to relate these turbulent correlations to the mean flow and thereby obtain a closed 
system. In the present analysis, the apparent mass flux term p'v', the apparent shea r 
stress term, pu’v’ (Reynolds stress term), and the apparent heat flux term Cppv’T' 
are modeled or represented by a new velocity component v, an eddy viscosity e , and an 
eddy conductivity Kq>, respectively. 

A new velocity component normal to the surface is defined as follows: 


v = 


v + 


p’v’ 


The eddy viscosity is defined as 


( 7 ) 


e = 


u'v’ 
p 9u/9y 


( 8 ) 
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and the eddy conductivity as 


v’T* 


( 9 ) 


The static turbulent Prandtl number is defined as follows: 


N 


_ u’v’ /9T/Sy\ 

Pr ’ t ~r¥\ */*%) 


Equation (10) can then be expressed in terms of equations (8) and (9) as 


( 10 ) 


Cp6 


(ID 


In terms of equations (7) to (11), the governing differential equations (eqs. (4), (5), 
and (6)) may be written as follows: 

Continuity: 


^-( r ip u ) + ^( r ipv) = o 


( 12 ) 


Momentum: 


pju 4^- + v 

ax 


^ = + 4 9 / rj . au\ 

ay/ dx r } d y\ 9 y/ 


(13) 


Energy: 


(“l;( c p T ) tf i( c p T ! 


dp -/3u\ 2 1 a 

= u — + ue — + -r 


dx 


\ay ) r 3 


rJjj. 


Np r 9y 


(CpT) 


(14) 


The terms e and e appearing in equations (13) and (14) are respectively defined as 
follows: 


= ( 1 + l r ) 


and 


1+ ±^EJL T] 


V Np r ,t 


(15) 


(16) 


The function r appearing in equations (15) and (16) represents the stream wise 
intermittency distribution in the transitional region of the boundary layer; T is a func- 
tion only of the x-coordinate and is discussed in detail subsequently. 

In order to complete the system of equations, the perfect gas law and Sutherland's 
viscosity relation are introduced: 
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Gas law: 


y - 1 

p = c p — — p t 


(17) 


Viscosity law: 

JL = (T_\ s/2 T e + s 
M e [r e ) T + S 


(Air only) 


(18) 


It should be noted that any viscosity relationship can be directly incorporated into the 
solution technique and that the system is not restricted to perfect-air (y = 1.4) boundary- 
layer flows. 


The system of governing equations then consists of three nonlinear partial differen- 
tial equations and two algebraic relations. Two of these differential equations (eqs. (13) 
and (14)) are second order and the remaining differential equation (eq. (12)) is first order. 
Consequently, if suitable relations for e, Np r and T can be specified, there are 
five unknowns, namely, u, V, p, T, and p and five equations. 


The pressure-gradient term appearing in equations (13) and (14) can be replaced 
by the Bernoulli relation; namely, 
dp du g 

dx - Pe u e ^ 


(19) 


which is determined from an inviscid solution. If variable entropy is considered, dp/dx 
is retained in equations (13) and (14). (See ref. 23.) 


Boun dary conditions .- In order to obtain a unique solution to the system of governing 
equations, it is necessary to satisfy the particular boundary conditions of the problem 
under consideration. These conditions are shown schematically in figure 2. ^Note that 
the outer edge conditions u e (x) and T e (x) are related through the energy equation 
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T|. e = T e + u e ^/ 2Cp^), where T t e is constant. The velocity and temperature distribu- 
tion at the edge of the boundary layer are determined from the shape of the body by using 
inviscid flow theory. The no-slip condition is imposed at the wall; however, arbitrary 
distributions of v w and T w or q w may be specified. 

The parabolic nature of equations (13) and (14) requires that the initial velocity and 
temperature profiles be specified at x^. These initial profiles are obtained in the pres- 
ent investigation from an exact numerical solution of the similar boundary-layer equa- 
tions for laminar flow (eqs. (B47) to (B49)). 


Transformed Plane 

The system of governing equations is singular at x = 0. The Probstein-Elliott 
(ref. 34) and Levy-Lees (ref. 35) transformation is utilized to remove this singularity as 
well as to reduce the growth of the boundary layer as the solution proceeds downstream. 
This transformation can be written as follows: 


£(x) = J PgU e ju e r^dx ( 20 a) 

r?(x,y) = -£• dy (20b) 

J 0 p e 

where the parameter t appearing in equation ( 20 b) is the transverse curvature term 
and is defined as 

t = 1 + ~ ( 21 a) 

r o 

or in terms of the y-coordinate as 


t = 


1 


+ 


COS <p 
o 


( 21 b) 


The relation between derivatives in the physical (x,y) and transformed (£, 77 ) coordinate 
systems is as follows: 



Pe u eMe r o ] 



> 

(a\ _ P e u e r o t3 /p\/ 8\ 

\ 9 y/x \j 2 l [Pel \ | 


( 22 ) 
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Two new parameters F and © are introduced and are defined as 


© = JL 


as well as a transformed normal velocity 

V = ^ (f + PVrpt (24) 

The governing equations in the transformed plane can then be expressed as follows: 
Continuity: 

lf +2 *f +F = ° (25) 


Momentum: 


Energy: 


2|F f + V |E - f (t% |£] + jj(f 2 - e) = o 

9£ d V dr l\ 9r? / 


2£F-^ + v|2- |-(t 2 > -L_ ? |2\ - oat 2 ief^ =0 

977 drj\ N Pr dri I UrJ 


where 


c p T e 

. 2£ du e 
u e d£ 


The parameter l can be written by using the viscosity relation (eq. (18)) and the equa- 
tion of state (eq. (17)) as 

1 - Kit!) < 29 


(for air only) where S = S/T e . 
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The transverse curvature term can be expressed in terms of the transformed 
variables as 


t 


± (l , 2 S cos ± 
\ p e u e 



(30) 


The physical coordinate normal to the wall is obtained from the inverse transformation; 
namely, 


y = 


cos <p 


-1 ± 



2 }/2£ cos cp 
Pe u e r o j 



(31) 


The positive sign is used in equations (30) and (31) for axisymmetric flow over bodies of 
revolution and the negative sign is used for flows inside axisymmetric ducts (nozzles). 

The boundary conditions in the transformed plane are as follows: 


Wall boundary: 


FU,0) = 0 
VU,0) = V W U) 

®te,0) = e w U) 1 



Edge conditions: 

F(l.1e) = A 

e(«. 1 e) = v 


(32a) 


(32b) 


The boundary condition at the wall for the transformed V-component can be related to the 
physical plane as 


V w = 


\j2i 


i r ] 
L e o 


P w v w 

P e u e 


(33) 


where the no-slip constraint has been imposed on equation (24). Note that the apparent 
mass flux term appearing in equations (4) to (6) is zero at the wall and that v w = v w . 
Therefore, equation (33) can be expressed in terms of the mass flux at the wall as 


w 


\J2$ p w v w 


Pe r ; 


e o 


3 Pp u 


e e 


(34) 
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Turbulent Transport Models 


The turbulent boundary layer is treated as a composite layer consisting of an inner 
and outer region as shown schematically in figure 3. (See Bradshaw, ref. 36.) 



Figure 3 •- Two-layer turbulent boundary- layer model. 


Inner region model .- The eddy viscosity model used for the inner region is based 
on the mixing-length hypothesis as developed by Prandtl (ref. 37) in 1925. The eddy vis- 
cosity for this region referenced to the molecular viscosity may be expressed as follows: 


-2 

/_e\ _ PZ 9u 

\M/i M 9y 


(35) 


where Z, the mixing length, may be written as 


l = Kjy 


(36) 


The value of Kj has been obtained experimentally and has a value of approximately 0.4, 
the value which will be used in the present analysis. However, Van Driest (ref. 38) con- 
cluded from an analysis based upon experimental data and the second problem of Stokes 
(ref. 39) that the correct form for the mixing length in the viscous sublayer should be as 
follows: 


l 


Kjy 



(37) 


where the exponential term is due to the damping effect of the wall on the turbulent fluc- 
tuations. The parameter A is usually referred to as the damping constant. The 
exponential term approaches zero at the outer edge of the viscous sublayer so that the 
law- of-the -wall region equation, as expressed in equation (36), is valid. The damping 
constant A is a strong function of the wall boundary conditions and is defined as 
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A=26„fe 


- 1/2 


(38) 


Equation (38) was originally obtained for incompressible, zero-pressure-gradient, 
solid- wall flow; that is, p = Constant, dp/dx = 0, and v w = 0. The relation has, how- 
ever, been applied to compressible flows where p and v are evaluated locally across 
the sublayer. In the analysis the mixing length is defined as follows: 


l =Kjy 



(39) 


where v is the average value of the kinematic viscosity taken over the viscous sublayer; 


N s i 



"i- 


The density and viscosity appearing in equation (38) are evaluated locally in 


the present paper; however, these quantities may be evaluated at the local wall boundary 
conditions if desired. Recent work at the Langley Research Center (ref. 40) has indi- 
cated that for hypersonic flow, the density and viscosity should be evaluated at the local 
wall boundary conditions. The eddy viscosity for the inner region referenced to the 
molecular viscosity can then be written as follows: 



pK^y 2 




1 - exp 



(40) 


where A is defined in equation (38). Cebeci (ref. 41) has recently attempted to account 
for the effects of both pressure gradient and wall mass flux on the damping constant. 

The reader interested in a discussion of this work is referred to the papers by Cebeci 
(ref. 41) and Harris (ref. 23). 


Outer region model .- The eddy viscosity in the outer region is based upon the 
Clauser (ref. 42) model. The ratio of the eddy viscosity to the molecular viscosity in 
the outer region can be expressed as follows: 

($„ ■ K 2 IT C <«> 


where 5* nc is the incompressible displacement thickness; 

6 inc= - F ) d y (42) 

J 0 

The use of 5* nc as the scaling parameter for the mixing length is discussed by Maise 
and McDonald (ref. 43). (See also Morkovin (ref. 44).) The value of K 2 in equation (41) 
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is taken to be 0.0168 as reported in reference 45. However, in order to account for the 
intermittent character of the outer layer flow, equation (41) is modified by an intermit - 
tency factor obtained by Klebanoff (ref. 46); that is, 


( e\ „ P u e ** - 

(m)o 2 M 6 inc y 

where the transverse intermittency factor y(y) is defined as 


(43) 


y = 


1 - erf 



2 


(44) 


The boundary -layer thickness 6 appearing in equation (44) is defined as the distance 
normal to the wall boundary where F = 0.995. 

Match ing procedure .- The criteria used to determine the boundary between the inner 
and outer regions is the continuity of eddy viscosity. A sketch of a typical eddy viscosity 
distribution is presented in figure 4. 



Figure 4.- Matching procedure for two-layer model. 

The matching procedure may then be formally written as follows: 



The location of the boundary separating the two regions y m is determined from the con- 
tinuity of equations (45); that is, where 
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Eddy conductivity .- The eddy conductivity is formulated in terms of a static turbu- 
lent Prandtl number Np r t and the eddy viscosity e. (See eqs. (9) to (11).) The two- 
layer concept for the eddy viscosity model suggests that there should also be a two-layer 
model for the static turbulent Prandtl number. Numerous assumptions have been made 
over the past years concerning the eddy conductivity, and corresponding models have been 
proposed to attempt to predict the mean turbulent boundary-layer temperature profiles. 
One of the earliest assumptions that has been used extensively is that the static turbulent 
Prandtl number is a constant equal to unity and implies that the heat and momentum are 
transferred by the same process. However, the data which are available, although often 
inclusive, definitely show that Np r> t is a function of y/5 for both incompressible and 
compressible flows. Consequently, the assumption of Np r> t = 1 would be expected to 
lead to error in predicting the temperature profiles for the turbulent boundary layer. 

The incompressible data which are available for the outer region of pipe flow 
(ref. 47) and boundary layers (ref. 48) indicate that Np r ^ ranges between 0.7 and 0.9. 
These data indicate that as the wall is approached, Np r t achieves a maximum value 
near the wall and then decreases rapidly to a value between 0.5 to 0.7 at the wall. (See 
also ref. 49.) Simpson, Whitten, and Moffat (ref. 50) found that Np rt ranged from 
approximately 0.95 at (y/6) » 0.1 to 0.45 at (y/5) = 1.0. The data in this region were 
predicted well by the expression Np r ^ = 0.95[l - 0.5(y/5)2j as proposed by Rotta 
(ref. 51). Cebeci (ref. 52) presents a continuous expression for the eddy conductivity 
which accounts for pressure gradient and mass transfer at the wall boundary. The 
resulting Np r t distribution appears to agree with both the experimental data of 
Simpson, Whitten, and Moffat (ref. 50) as well as with the Jenkins model (ref. 53). For 
compressible flows, where very little data are available, it appears that Np r ^ has a 
value near unity in the outer region of the boundary layer and a value between 0.7 and 0.9 
at the wall boundary (ref. 49). Rotta (ref. 51) found that Np r ^ may reach values as 
high as 2.0 as the wall is approached before decreasing rapidly to the wall boundary value. 
Meier and Rotta (ref. 54) found for 1.75 £ Moo = 4.5 that Np r t increased above unity 
for y + < 50 and ranged between 0.8 and 0.85 as the outer edge of the boundary layer was 
approached. 

The eddy conductivity and resulting static turbulent Prandtl number model must of 
necessity at the present time be formulated in terms of existing experimental data. The 
assessment of the value of the current empirical models must be based upon the agree- 
ment between experimental and calculated temperature profiles over a wide range of flow 
variables. The inconclusiveness of much of the experimental data which exist to date for 
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both subsonic and supersonic flows and the lack of data for hypersonic flows clearly 
indicates the severe need for accurate, well-documented experimental temperature and 
velocity profile data. In the present analysis, since the cases considered deal specifically 
with either wall gradient values or velocity profile parameters, a constant Np r ^ equal 
to 0.9 is utilized. However, the numerical procedure utilized is completely general and 
any function Np r ^(y/S) could be used. 

Reynolds number effects .- The eddy viscosity model as presented in the previous 
sections has implied that Kj, K 2 , and A + (A + = Au T ,wAw) are constants (see 
eqs. (45)) and therefore independent of flow conditions. However, recent studies have 
shown that Kq, K 2 , and A + are strong functions of R e> £ for moderate to low 
Reynolds numbers; for example, R g q < 5000. (See refs. 52 and 55.) Cebeci (ref. 52) 
has found that a Kq value of 0.4, as used in the present analysis, is accurate and reli- 
able in predicting the velocity profiles in the law- of -the -wall region provided that 
R e ,e > 6000. Cebeci, using the data of Simpson (ref. 56), generated polynomial curves 
for Kq and A + as functions of R g q. These functions were then used to obtain pre- 
dictions for the data of Whitten (ref. 57) which were in very good agreement in the law- 
of-the-wall region for 1200 s R e ^ ^ 4100. McDonald (ref. 55) has shown that K 2 is a 
strong function of Reynolds number for R e q < 5000; that is, K 2 ranges from a value 
near that used in the present analysis for R e q > 5000 to values as large as 0.036 for 
^e,0 = ^50. 

For all the cases presented in the present paper, Kq, K£, and A + are held 
constant and equal to 0.4, 0.0168, and 26, respectively, since low R e q cases are not 
considered. However, the flexibility of the numerical technique allows any desired varia- 
tion of these variables to be easily incorporated into the solution. 

Transformed models .- Since the governing equations are solved in the transformed 
plane, it is necessary to transform the eddy viscosity relations from the real plane to the 
transformed plane. 


In the inner region the ratio of eddy viscosity to molecular viscosity is as follows: 
I e\ p e u e^l r o + j y^n^ 

r V 


9F 


d V 


(47) 


where y is defined by equation (31). The parameter n appearing in equation (47) is 
the damping term and is defined as 


n = 1 - exp^njllg) 

where 


(48) 


n l 


0 » l 


(49) 
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and 


U 2 = 


yp* u 


e e 


26Z0 2 ju 
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^e r c/w / 8F 


eL 


VZ? Ww 


1/2 


(50) 


In the outer region the ratio of eddy viscosity to molecular viscosity is as follows: 


where 


(&- 


p e Ue v r6 inc 
K 2 — 




Z0 z 


r = 


1 - erf 5( — - 0.78 


(51) 


(52) 


and 


5 inc = -®r I ^ “ F ) 


P e u e r o 0 


d?7 


(53) 


Transition Region 

Equations (25), (26), (27), and (29) together with the boundary conditions (eqs. (32)), 
and the eddy viscosity relations defined by equations (47) and (51) complete the required 
system for either laminar or fully developed turbulent boundary -layer flows. However, 
the main objective is to present a technique that will efficiently solve the laminar, transi- 
tional, or turbulent boundary-layer equations as the boundary layer develops along the 
surface. Consequently, the location of transition Xo. ,, the extent of the transitional flow 

l 

x,. f - x, . , and the characteristics of the mean flow structure in the transition region 

i,i t,i 

must be taken into consideration. 

Stability and transition .- Stability theory cannot currently be used to predict either 
the nonlinear details of the transition process after the two-dimensional waves have been 
amplified or the location of transition x^. Stability theory can, however, establish the 
unstable boundary-layer profiles and the initial amplification rates. The theory can 
identify those frequencies which will be amplified at the greatest rate as well as the effect 
on stability of various flow parameters. (See refs. 2, 58, and 59.) A review of methods 
used to predict the location of transition from stability theory is presented by Jaffe, 
Okamura, and Smith in reference 60. However, the connection, if any, between stability 
and transition has not yet been established. 

Transition location .- Many parameters influence the location of transition. These 
parameters can best be thought of as forming a parameter phase space. Such a param- 
eter phase space would include Reynolds number, Mach number, unit Reynolds number, 


24 



surface roughness, nose bluntness, pressure gradients, boundary conditions at the wall, 
angle of attack, free-stream turbulence level, and radiated aerodynamic noise. Morkovin 
(ref. 12) recently completed a very extensive and thorough examination of the current 
state of the art of transition from laminar to turbulent shear layers. The most striking 
conclusion that one obtains from the review is that although a great bulk of experimental 
data on transition currently exists, much of the information on high-speed transition has 
not been documented in sufficient detail to allow the separation of the effects of these 
multiple parameters on transition. A discussion of the effects of each of the param- 
eters that may influence transition in high-speed flow is beyond the scope of the present 
paper. 

The reader interested in a detailed discussion is directed to the paper by Morkovin 
(ref. 12) where over 300 related references are cited and discussed. Another, although 
less detailed discussion, is presented by Fischer (ref. 61). The effects of radiated aero- 
dynamic noise on transition are discussed by Pate and Schueler (ref. 62). Hypersonic 
transition, to name but a few references, is discussed by Softley, Graber, and Zempel 
(ref. 63), Richards (ref. 64), Potter and Whitfield (ref. 65), and Deem and Murphy 
(ref. 66). A discussion on the effects of extreme surface cooling on hypersonic flat- 
plate transition is presented by Cary (ref. 67). The effects of nose bluntness and surface 
roughness on boundary-layer transition are discussed by Potter and Whitfield (ref. 68). 

In the present analysis the location of transition will be determined by one of, or a com- 
bination of, the following three methods. These methods are (1) use of a stability index 
or vorticity Reynolds number first proposed by Rouse (ref. 69), (2) use of correlations 
based upon a collection of experimental data over a broad range of test conditions, and 
(3) use of the measured experimental location of transition as a direct input into the 
analytical solution. 

Stability index .- Hunter Rouse (ref. 69) nearly 25 years ago obtained by dimensional 
analysis a stability index expressed as follows: 


_ y 2 8u 
x " “8y 


(54) 


This index has the form of a vorticity Reynolds number which is obtained from the ratio 
of the local inertial stress py2(9u/8y)^ to the local viscous stress p(9 u/8y). Rouse 
assumed that for transition to occur the stability index should reach some limiting value 
which was assumed to be invariant. He was able to show further that this invariant value 
(X ma x) cr s h° u ld on the order of 500 for incompressible flows. 

The use of Xmax as a stability index is in principle similar to the basic reasoning 
which led Osborne Reynolds (ref. 31) to postulate that the nondimensional parameter ud/V 
could be used to define a critical value (ud/i/) cr at which transition would occur in a 
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circular pipe of diameter d. Unfortunately, (Xmax)cr is a ^ unction °f transition 
parameter phase space in much the same fashion as the critical Reynolds number and 
cannot in reality be a true invariant for transition as suggested by Rouse. The stability 
index does, however, possess a number of important characteristics which can be 
directly related to the stability of laminar flows. 

A schematic distribution of x for a compressible laminar boundary layer is pre- 
sented in figure 5. The vorticity Reynolds number has a value of zero at the wall and 
approaches zero as the outer edge of the layer is approached. The maximum value of 
X, that is, Xmax’ occurs at some transverse location (y/ 5 )x max ‘ The values of X max 

and (y/5) v are of importance in the use of the vorticity Reynolds number as a 
X m ax 

guide to boundary-layer stability and transition. 

As the laminar boundary layer develops over a surface, Xmax increases mono- 
tonically until the critical value (Xmax)cr is reached at which point transition is 
assumed to occur; that is, the location of x^. For compressible flows (Xmax) cr is 
not an invariant. In the present study (x max ) cr ranged from approximately 2100 to 
values on the order of 4000. The variation (x ma x)cr is a strong func t ion of unit 
Reynolds number for data obtained in air wind-tunnel facilities. However, although not 
invariant, the stability index does exhibit the same dependence on various parameters 
that is predicted by the more complicated stability theory. For example, at a given 
streamwise location x, the value of x m ax for su P ersonic flow is found to decrease 
(which implies a more stable flow) with wall cooling, wall suction, and favorable pres- 
sure gradients, whereas it increases (which implies a more unstable flow) with wall 
heating, mass injection (transpiration), and adverse pressure gradients. The vorticity 
Reynolds number X max appears to have been used as a correlation parameter in Only 
two boundary- layer transition studies. A modified form of the parameter was used to 



Figure 5 »- Vorticity Reynolds number. 
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correlate the effect of free -stream turbulence on transition by Van Driest and Blumer 
(ref. 70). Correlation attempts, using Rouse's original invariant assumptions, were made 
by Ivensen and Hsu (ref. 71); however, the results were only fair. 

One of the most important characteristics of the vorticity Reynolds number is that 
the value of (y/^)x max * s exce H en t agreement with the experimental location of the 
critical layer which represents the distance normal to the wall where laminar flow break- 
down will occur. Stainback (ref. 72) computed the Rouse stability index for similar 
laminar-boundary-layer flows over a broad range of ratios of wall temperature to total 
temperature for Mach numbers up to 16. The numerical calculations were made for both 

air and helium boundary layers. The agreement between (y/6) v and the experi- 

*max 

mental critical layer position was excellent over the entire range. 

Empi r ical correlations .- In most instances one has to rely on empirical correla- 
:ions of experimental transition data in order to fix the most probable transition location 
’or a design configuration. However, caution should always be used when obtaining the 
nost probable location of transition from such correlations, since any given correlation 
s based upon a specific collection of data which will not be completely general. There 
:urrently exists a number of empirical correlations for predicting the probable location 
>f transition. Some of these correlations are of questionable value; however, some can 
ie used with confidence providing it is realized that one is predicting a probable range 
if locations and not an exact fixed point. One of the more successful correlations 
unpublished) was obtained by Stainback and Beckwith for sharp cones at zero angle of 
ittack. The correlation is based on experimental transition data obtained over a wide 
range of test conditions in air wind tunnels, ballistic ranges, and free flight. The correla- 
;ion can be expressed as follows: 


( /rs N 0.7exp(-0.05M|) 

lo gl0 —5— = p + QM e (© w ) V 67 


(55) 


where R„ j and R „* are the Reynolds numbers based on base diameter d and the 
C ^ e ? °t 

displacement thickness at transition, respectively. The constants P, Q, and R are 
functions of the environment in which transition was measured and are given in the fol- 
lowing table: 


Facility 

P 

Q 

R 

Air wind tunnel 

1.6217 

0.11708 

0.25 

Ballistic range 

2.2082 

.14047 

.15 

Free flight 

1.2683 

.10441 

.325 
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as 


Equation (55) can be expressed in terms of the transition Reynolds number 
follows: 



R c 


3 R e5 10 


2P + 2QM e (0 w ) 
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(0;094Mg + 1.220 w )' 


(56) 


Equations (55) and (56) are similar to the relations presented by Bertram and Beckwith 
(ref. 73); however, the constants P, Q, and R presented in the preceding table are 
somewhat different from those of reference 73. 


Experimental transition .- Much of the confusion that exists in the literature con- 
cerning boundary-layer transition may be attributed to the following two factors. The 
first factor is that in many instances the investigator who made the experimental study 
may not have carefully measured or recorded the complete conditions under which the 
data were obtained. The second factor is that the experimentally observed transition 
location depends on the experimental technique used. 


The importance of carefully measuring the environment under which the experi- 
ments are made cannot be overstressed. In the past, many of the factors which influence 
transition such as free-stream turbulence and acoustic radiation from the tunnel side- 
wall boundary layer were not measured. (See refs. 73 to 76.) The location of transition 
as obtained experimentally is a function of the method used to determine its location. 
There are currently a number of techniques used to obtain the transition location. Some 
of these methods are hot-wire traverses, pitot-tube surveys near the wall, visual indica- 
tion from schlieren photographs, and heat -transfer measurements at the wall. Each of 
these methods basically measures a different flow process. Consequently, it would be 
misleading to believe that each technique would yield the same location for transition if 
simultaneously applied to the same boundary layer. Of course, the concept of a transition 
"point" is misleading in itself since transition does not occur at a "point" but instead over 
some finite region. 


For the test cases presented herein, the transition location Xf. • will be determined 
from heat-transfer measurements at the wall whenever possible. The reason for this 
choice is that it is the method most often used in the literature. However, it should be 
noted that the actual nonlinear transition process begins somewhat upstream of the loca- 
tion where the heat transfer at the wall deviates from the laminar trend. 


Transitional flow structure .- Once the transition location has been fixed for a given 
problem, one must next consider the following two important factors; first, the length of 
the transition region, ^ f ~ x t j (sometimes referred to as transition extent), and sec- 
ondly, the mean flow characteristics within the region. Once appropriate models are 
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obtained for these two factors, it is possible to connect the three flow regions smoothly 
so that one set of governing equations may be used. 

The laminar boundary -layer equations (eqs. (1) to (3)) should yield accurate profiles 
and wall boundary gradients in the region prior to turbulent spot formation. The inter- 
mittent appearance of the turbulent spots and the process of cross-contamination are not 
well understood. The spots originate in a more or less random fashion and merge with 
one another as they grow and move downstream. Eventually, the entire layer is con- 
taminated which marks the end of the transition process Xj. As the turbulent spots 
move over a fixed point in the transition region, the point experiences an alternation 
between fully laminar flow when no spot is present and turbulent flow when engulfed by a 
spot. These alternations can be described by an intermittency factor which represents 
the fraction of time that any point in the transition region is engulfed by turbulent flow. 


The distribution of spots in time and space is Gaussian for low-speed natural 
transition. However, very little is known about the spot distribution in high-speed com- 
pressible flow. Furthermore, there is no assurance that the spot formation and distribu- 
tion in hypersonic flows will be analogous to the low-speed model; however, in the absence 
of a more satisfactory theory, the approach of Dhawan and Narasimha (ref. 77) is used. In 
reference 77 the source density function of Emmons (ref. 20) was used to formulate the 
probability distribution (intermittency) of the turbulent spots as follows: 


F(i) = 1 - exp 



(57) 


where 


1 = 



(58) 


for x^. ^ § x = Xj j. The term f in equation (57) represents a normalized streamwise 
coordinate in the transition region, and X is a measure of the extent of the transition 
region; that is, 

X = ( x )p_3 / /4 “ ( x )p = i/4 (59) 


The intermittency distributions in the transverse direction (y-direction) are similar to 
those observed by Corrsin and Kistler (ref. 78) for fully developed turbulent boundary 
layers (ref. 79). Corrsin and Kistler found that the transverse intermittency varied 
from a maximum of unity near the wall to a near-zero value at the outer edge of the layer. 
The transverse intermittency distribution is of secondary importance in relation to the 
streamwise distribution in determining the mean profiles and wall fluxes; consequently, 
the only intermittency distribution applied in the transverse direction (y-direction) is that 
of Klebanoff (ref. 46) in the outer layer as applied to the fully developed turbulent layer. 
(See eq. (44).) It should be carefully noted in order to avoid any confusion concerning the 
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two intermittency functions used herein that r is a function only of the x-coordinate 
whereas y is a function only of the y-coordinate. 

Transition extent .- The assumption of a universal intermittency distribution implies 
that the transition zone length (transition extent) can be expressed as a function of the 
transition Reynolds number u e x£ ^jv e . In reference 77 it is shown, for the transition 
data considered, that the data are represented on the average by the equation 


R, 


! > Ax t " 5 ( Re » x t,i) 


0.8 


(60) 


where R e ,Ax t = Ax t- 


The location of the end of transition 


directly from equation (60) as follows: 


*t,f 


can then be obtained 


= + **?{***$* (61 > 
where R e is the local unit Reynolds number u e /^ e . Morkovin (ref. 12) noted that only 
about 50 percent of the experimental data he considered could be fitted to the low-speed 
universal curve of Dhawan and Narasimha, that is, to equation (60). This result was to 
be expected, since the data considered in reference 77 covered only a very limited Mach 
number range. 


Potter and Whitfield (ref. 68) measured the extent of the transition zone over a 
rather broad Mach number range (3 t Moo § 5; M w = 8). They observed that the transi- 
tion region, when defined in terms of R e ^ x ^_, is basically independent of the unit Reynolds 
number and leading-edge geometry; that is, 

^e,AXj. = ^^e,Xj. ^ ( 62 ) 

They noted (ref. 68) that the extent of the transition region increased with increasing 
transition Reynolds number over the Mach number range 0 £ M^ = 8 for adiabatic walls. 
The extent of the transition region was also observed to increase with increasing Mach 
numbers for a fixed transition Reynolds number. 


In the present analysis, because of the lack of general correlations for the extent of 
transition, this quantity is obtained directly from the experimental data unless otherwise 
noted. In particular, if heat -transfer data are available, the transition zone is assumed 
to lie between the initial deviation from the laminar heat -transfer distribution and the 
final peak heating location. The transition region defined on the basis of the Stanton num- 
ber distribution is presented in figure 6. The design engineer seldom has the advantage 
of having access to experimental data which were obtained under the actual flight condi- 
tions. Consequently, the most probable location of transition would be obtained from a 
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Figure 6.- Transition extent definition. 


correlation such as presented in equation (56). The extent of transition could then be 
obtained from equation (61) or an approximate relation such as follows (see ref. 80): 


Re »*t.f _ „ 

R e x t ■ 
e > x t,i 


(63) 


NUMERICAL SOLUTION OF THE GOVERNING EQUATIONS 

The system of governing equations for compressible laminar, transitional, and 
turbulent boundary layers consist of three nonlinear partial differential equations (see 
eqs. (25) to (27)) and two algebraic relations. (See eqs. (17) and (18).) The most impor- 
tant feature of this system is that it is parabolic and therefore can be numerically inte- 
grated in a step-by-step procedure in the streamwise direction. In order to cast the 
equations into a form in which the step-by-step procedure can be efficiently utilized, the 
derivatives with respect to £ and 77 are replaced by linear finite -difference quotients. 

The method of linearization and solution used in the present analysis closely paral- 
lels that of Fliigge-Lotz and Blottner (ref. 29) with modifications suggested by Davis and 
Fliigge-Lotz (ref. 30) to improve the accuracy. These modifications involve the use of 
three-point implicit differences in the ^-direction which produce truncation errors of 
order A£j rather than Ax as in reference 29. The primary difference between 

the present development and that of reference 30 is that the solution is obtained in the 
transformed plane for arbitrary grid-point spacing in the ^-direction and for a spacing in 
the 77 -direction such that the ratio of the spacing between any two successive grid points 
is a constant. 
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The Implicit Solution Technique 


Finite -difference mesh model .- It has been shown for laminar boundary layers that 
equally spaced grid points can be used in the normal coordinate direction. (See refs. 29 
and 30.) However, for transitional and turbulent boundary layers, the use of equally 
spaced grid points is not practical because the fine-mesh size required to obtain accurate 
results near the wall boundary is inefficient for the entire boundary layer. The grid- 
point spacing in the 77 -direction used in the present analysis assumes that the ratio of any 
two successive steps is a constant; that is, the successive A 77 ^ form a geometric pro- 
gression. (See ref. 81.) 

The advantage of having variable grid-point spacing in the x-coordinate becomes 
clearly apparent for problems in which either the rate of change of the boundary condi- 
tions is large or discontinuous or the mean profiles are changing rapidly. A good example 
would be a slender blunted cone in supersonic flow. Variable grid-point spacing in the 
^-direction would be utilized by having very small steps in the stagnation region where 
the pressure gradient is severe (favorable) and again in some downstream region where 
transitional flow exists. A good example of discontinuous boundary conditions would be 
a sharp-tipped cone with a porous insert at some downstream station through which a gas 
is being injected into the boundary layer. Relatively large step sizes could be utilized 
upstream of the ramp injection; however, small steps must be used in the region of the 
ramp injection. Downstream of the porous region, as the flow relaxes, larger step sizes 
could be used. It is very important that small grid-point spacing be utilized in the transi- 
tion region where the mean profiles are a strong function of the intermittency distribution. 

In constructing the difference quotients, the sketch of the grid-point distribution pre- 
sented in figure 7 is useful for reference. The dependent variables F, 0, and V are 
assumed to be known at each of the N grid points along the m-1 and m stations, 
but unknown at station m + 1. The A£j and A £ 2 values, not specified to be equal, 
are obtained from the specified x-values x m , and and equation ( 20 a). 

The relationship between the A 77 ^ for the chosen grid-point spacing is given by the fol- 
lowing equation: 

A7 ?i = (K ) i_1 At ?1 (i = 1,2,. . ., N) (64) 

where K is the ratio of any two successive steps, A? 7 j is the spacing between the 
second grid point and the wall (note that the first grid point is at the wall boundary), and 
N denotes the total number of grid points across the chosen 77 -strip. The total thick- 
ness of the 77 -strip can then be expressed as follows: 


% 


= A77j 


1 - (k)N- 1 

1 - K 


(K * 1) (65) 
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Figure 7*- Finite-difference grid model. 


The selection of the optimum K and N values for a specified depends upon the 
particular problem under consideration. The main objective in the selection is to obtain 
the minimum number of grid points with which a convergent solution may be obtained and 
thereby minimize the computer processing time per test case. The laminar boundary 
layer presents no problem since a K value of unity is acceptable; however, for transi- 
tional and turbulent layers, the value of K will be a number slightly greater than unity, 
for example, 1.02 to 1.04. If transitional or turbulent flow occurs in a given problem; 
the laminar part of the boundary layer is calculated with the value of K used for the 
turbulent region; that is, for a given problem K is invariant. 

Difference equations .- Three-point implicit difference relations (see appendix A) 
are used to reduce the transformed momentum and energy equations (eqs. (26) and (27)) 
to finite -difference form. The difference quotients produce linear difference equations 
when substituted into the momentum and energy equations provided that truncation terms 
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of the order A| m _i A£ m and A?7 n _ 1 A % are neglected. The resulting difference 
equations may be written as follows: 


A (1) -o' , td( 1) ■p' . /'•(l)-p’ . T)^)© 

^n F m+l,n-l + n F m+l,n + C n F m+l,n+l + D n °m+l,n-l 


+ E i 1)0 m + l,n + 4 1)e m + l,„ + l = 4 1) 


( 66 ) 


■^n ^ F m+l,n-l 


B^F * + C^F , 1 

n m+1 ,n n m+l,n+l n 


J m+1 ,n-l 


E^0 < + F^0 . . = G^ 

n m+l,n n m+l,n+l n 


(67) 


The coefficients . . ,,G^, A^, . . ,,G^ (see appendix B) are functions 

of known quantities at stations m and m - 1. The dependent variable V does not 
appear explicitly as an unknown at station m + 1 in these equations (eqs. ( 66 ) and (67)). 
This condition occurs because of the particular way that the nonlinear terms V(3F/ 3 ij) 
and V( 30/ 877 ) appearing in the momentum and energy equations (eqs. (26) and (27)), 
respectively, are linearized (see eq. (A23)); however, equations ( 66 ) and (67) are coupled 
in that the dependent variables F and 0 appear explicitly in each equation. 

Solution of difference equations .- The system of difference equations (eqs. (66) 
and (67)) represents a set of exactly 2(N - 1) equations for 2(N - 1) unknowns since 
the N-l unknown values of V do not appear explicitly and the A^,B^, . . ,,G^ 
coefficients are known. The proper boundary conditions to be used with the difference 
equations are specified in equations (32). The system of algebraic linear equations may 
be written as a tridiagonal matrix; consequently, an efficient algorithm is available for 
their simultaneous solution. The simultaneous or coupled-solution technique has been 
discussed in detail by Flugge-Lotz and Blottner (ref. 29). For perfect gas flows or sys- 
tems where only two governing equations must be considered, the simultaneous solution 
technique (coupled solution) is very efficient since the equations are solved without the 
need for iteration procedures such as in reference 33. However, for real-gas flows or 
systems where more than two governing equations are required, the iterative (uncoupled) 
solution technique proposed by Blottner (ref. 6) should be used to avoid the large digital 
computer storage requirements associated with the present coupled solution technique. 
Another uncoupled solution technique that appears to be more flexible and faster than the 
method presented in reference 6 is that developed by Keller and Cebeci (ref. 82). This 
method should efficiently handle both three-dimensional steady and two-dimensional 
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unsteady flows as well as flows with chemically reacting species; however, the technique 
to date has not been used extensively. 

Solution of continuity equation. - The continuity equation (eq. (25)) can be numerically 
solved for the N - 1 unknown values of V at station m + 1 once the values of F 
and 0 are known at station m + 1. Equation (25) is integrated to yield the following 
relation for V at the grid point (m + 1, n): 


^m+l,n ^rn+1,1 



( 68 ) 


where V m + l,l represents the boundary condition at the wall and is defined in equa- 
tion (34) as a function of the mass transfer at the wall (pv) w . The integral appearing in 
equation (68) is numerically integrated across the 77-strip to obtain the N - 1 values of 
V. In the present analysis, the trapezoidal rule of integration was utilized; however, any 
sufficiently accurate numerical procedure could be used. 


Initial profiles .- Initial profiles for starting the finite-difference scheme are 
required at two stations since three-point differences are utilized. The initial profiles 
at the stagnation point or line for blunt bodies, or near x = 0 for sharp-tipped bodies, 
are obtained by a numerical solution of the similar boundary-layer equations. The equa- 
tions are solved by a fourth-order Runge-Kutta scheme with a Newton's iteration method 
to modify the initial estimates of the gradients of F and © evaluated at the wall bound- 
ary. The N-l values of F, 0, and V which are now known at the N - 1 equally 
spaced grid points are numerically redistributed to N-l grid points whose spacing is 
determined from equations (64) and (65), if a variable spacing is required. The second 
initial profile located at station m is assumed to be identical to the one located at sta- 
tion m - 1. Any errors that might be incurred because of this assumption are minimized 
by using an extremely small A£, that is, an initial step size in the physical plane on the 
order of Ax = 1 X 10“ 5. The solution at the unknown station m + 1 is then obtained by 
the finite-difference method. Extremely small, equally spaced A£-steps are used in the 
region of the initial profile. The step size is increased after the errors due to the starting 
procedure have approached zero, that is, after 10 to 15 steps in A£. 

It is also advantageous to have the capability of starting the solution from experi- 
mentally measured profiles, especially in the case of turbulent flow. This capability 
can be directly incorporated into the numerical procedure used in the present analysis. 

This capability is extremely useful for cases where one cannot easily locate the origin 
of the boundary layer, for example, on nozzle walls. 


Evaluati on of wall derivatives .- The shear stress and heat transfer at the wall are 

F and 0 evaluated at the wall, respectively. 


directly proportional to the gradient of 
By using G to represent a general quantity where 


> m+ i 1 is n °t specified to be zero, 
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the four-point difference scheme used to evaluate derivatives at the wall is as follows 
(see pp. 49 to 50, ref. 29): 




Vm+l , 1 _ Y 7 G m+l,l + Y 8 G m+l ,2 + Y 9 G m+l,3 + Y 10 G m+l,4 
where the coefficients Y^, . . Y^g are defined by the following relations: 


= (l + K + K 2 ) 2 [jC(l +K)-l| + l+ K 
7 (1 + K)(l + K + K 2 )K 3 A?7 1 




y 8 = 


Y 9 = - 


Y 10 = 


1 + K + Kf 

K 2 A77j 

1 + K + K 2 
(1 + K)K 

1 

(l + K + K 2 )K 3 A?7 1 


( 69 ) 


(70) 




For the case of equally spaced grid points in the ^-direction (K = 1) , equations (70) 
become 


Y 7 = - 

' 6 A 77 


*N 


Y 0 = 


18 


8 6 A?? 


Y = - g 
9 6 A 77 


(71) 


Y i0 = 


10 6 Arj 




and equation (69) reduces to the four-point relation (see pp. 49 to 50, ref. 29); that is, 

(f) m+M - - 6 ^( 11 G m + l,l ' 18 <W,2 + 9 G m t l,3 ' 2 G m +M ) < 72 > 
Eddy Viscosity Distribution 

The primary difference between the present solution technique and other typical 
procedures such as those of Cebeci, Smith, and Mosinkis (ref. 83) and Beckwith and 
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Bushnell (ref. 84) is that the momentum and energy equations (eqs. (66) and (67)) are 
simultaneously solved without iteration, whereas in the techniques of these two refer- 
ences the momentum and energy equations are each individually solved and iterated for 
convergence. In the present procedure, extreme care must be used when the eddy vis- 
cosity functions e and e (see eqs. (15) and (16)) and their derivatives with respect to 
77 are extrapolated from the known values of e m _^ n and i m ,n to the unknown sta- 
tion m + 1, n. 

During the development of the present digital computer program, the calculations 
would frequently become unstable in either the transitional or turbulent flow region. 

This problem would always occur in one of two ways. In some instances an apparently 
converged solution would be obtained, but the distribution of boundary- layer thickness 
would not be smooth. In other instances, where the transition was abrupt or where bound- 
ary conditions were abruptly changed, the solution would not converge. The problem was 
caused by a wavy or rippled distribution in eddy viscosity across the layer. These rip- 
ples first occurred in the region where the inner and outer eddy viscosity models were 
matched. If the initial ripples were below a certain level, the solution would apparently 
converge, but slightly nonsmooth boundary-layer thickness distributions would occur. If 
the initial ripples were above a certain level, the oscillations would grow very rapidly 
as a function of $ and propagate throughout the layer as the solution proceeded down- 
stream with the result that valid solutions could not be obtained downstream of the initial 
ripples. 

The initial extrapolation of the known values of e m n and n to the unknown 

grid point (m + 1, n) is obtained as follows (see eq. (A3)): 

e m+l,n = X 4 e m,n " X 5 e m-l,n ( 73 ) 

However, there is no assurance that the distribution of the extrapolated values at sta- 
tion m + 1 will be smooth across the layer for all possible flow conditions. If oscilla- 
tions occur in the extrapolated e distribution of sufficient magnitude to cause the sign 
of the derivative of i with respect to 77 to alternate, then the entire calculation 
becomes highly unstable. 

The requirement of small grid-point spacing in the law-of -the -wall region con- 
tributes to the stability problem since the size of an "acceptable oscillation" is apparently 
a function of the grid-point spacing being utilized in the 77 -direction. For turbulent layers 
where the viscous sublayer is relatively thick, the grid-point spacing in the outer portion 
of the law-of-the-wall region will be large in comparison with cases where the viscous 
sublayer is relatively thin. Consequently, the former case can tolerate a larger ripple 
in the region of the match point y m than can the latter case without experiencing a 
change in the sign of the derivative of e with respect to 77 . 
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There are two possible ways to eliminate the problem caused by the oscillations 
in the eddy viscosity distribution. The first approach is to develop an iteration scheme; 
in which case the present solution technique has no advantage over other techniques (for 
example, refs. 33, 83, and 84), that is the advantages of the simultaneous solution would 
be lost. The second approach is to smooth numerically the extrapolated eddy viscosity 
distribution prior to the matrix solution. Both approaches were tried during the develop- 
ment phase of the digital computer program. The second approach was incorporated into 
the solution and is discussed in the remaining part of this section. 

The problem posed by the ripples in the eddy viscosity distribution, if they exist, 
can be avoided by utilizing a three-point mean value for e at station m + 1, n; that is, 

r- \ _ e m+l,n-l + e m+l,n + e m+l,n+l , nA \ 

' av 'm+l,n 3 K ’ 

where e av denotes the three -point mean of the eddy viscosity function. In the present 
procedure the eddy viscosity functions appearing on the right-hand side of equation (74) 
are first obtained at each grid point across the m + 1 station from equation (73). After 
these values are obtained, the three-point mean is evaluated at each of the N - 1 grid 
points from equation (74). The matrix solution for F and 9 is then obtained for 
equations (66) and (67) by using these precomputed values of (e av ) m+ i n . After the 
N - 1 values for F, 0, and V at station m + 1 have been obtained, the eddy vis- 
cosity distribution is recalculated at the m + 1 station from equations (47) and (51) 
prior to moving to the next £ grid-point station. This procedure has been found to be 
stable under all circumstances and to yield convergent solutions for transitional and fully 
turbulent boundary layers. 

EXAMPLE SOLUTIONS AND COMPARISONS WITH 
EXPERIMENTAL DATA 

The selection of a typical set of test cases always presents a problem since there 
are many possibilities from which to choose. However, the cases considered in the pres- 
ent paper were chosen to provide an indication of the merits of the solution technique as 
well as the validity of the eddy viscosity and intermittency models. In all cases pre- 
sented herein, the gas is assumed to be perfect air with a constant ratio of specific heats 
(y = 1.4), a constant Prandtl number (Np r = 0.72), and a constant static turbulent Prandtl 
number (Np r t = 0.9). The molecular viscosity n is evaluated from Sutherland's vis- 
cosity law (eq. (18)). The external pressure distributions used are either experimental 
or obtained from an exact solution of the full inviscid Euler equations. All calculations 
were made on the Control Data Corporation 6600 digital computer. 


38 



Flows with large adverse pressure gradients are not considered in the present 
paper. This omission is due to the current lack of understanding of the effect of large 
adverse pressure gradients on the turbulent models for supersonic flow as discussed by 
Harris (ref. 23) and Beckwith (ref. 17) who have made comparisons of numerical results 
with the experimental data of McLafferty and Barber (ref. 85). 

High Reynolds Number Turbulent Flow 

The accurate prediction of boundary-layer characteristics for high Reynolds num- 
ber turbulent flow is important in the design of high-speed vehicles. In particular, it is 
important to be able to predict with accuracy the skin-friction drag. A good example of 
high Reynolds number turbulent flow is provided by the data of Moore and Harkness 
(ref. 86). The experimental skin-friction data were measured with a floating-element- 
type balance. The data were obtained on a sharp-leading-edge flat-plate model that com- 
pletely spanned the width of the tunnel test section for 1 x 10^ < R e x < 2 x 10® and on 
the tunnel diffuser floor for 2 x 10® < Re,x < ^ x 10®* The test conditions were as 
follows: 

Mco = 2.8 

p t ^ = 0.997 MN/m 2 
T t,oo= 311.1 K 



The experimental transition location was not reported in reference 86. Conse- 
quently, for the numerical calculations, the transition location was determined by the 
stability index and was assumed to occur at the x-station where Xmax ac hieved a value 
of 2500. The extent of the transition region was calculated from equation (61). The 
intermittency distribution was calculated from equation (57). The solution was started 
at the leading edge of the sharp flat plate (x = 0) by obtaining a numerical solution of the 
similar boundary-layer equations. (See eqs. (B47) to (B49).) The grid-point spacing 
was varied in both the ij- and ^-directions in order to check for convergence. 

The numerical results for the skin-friction coefficient distribution are compared 
with the experimental data in figure 8(a). The agreement is excellent over the entire 
Reynolds number range of the experimental data for K = 1.04; however, for K = 1.01, 
convergence was not attained. It should be noted at this point that the terms convergence 
and stability, as used in relation to the numerical method, are defined herein as in refer- 
ences 23 and 29. The divergence for K = 1.01 is attributed to an insufficient number 
of grid points in the wall region, in particular, in the thin viscous sublayer region. 
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The effect of the grid-point spacing parameter K on the numerical solution was studied 
for additional values of 1.02, 1.03, 1.05, and 1.06. The solution was found to diverge for 
K < 1.02 and converge for K s 1.02. The C f e results for K & 1.02 (other than 
K = 1.04) are not presented since the variation of Cf e with K would not be discernible 
if plotted to the scale of figure 8(a). The convergence criteria used in this particular 
example was to decrease the rj-grid-point spacing until any change which occurred in C^ e 
at a given x-solution station was beyond the fourth significant digit. The laminar curve 
shown in figure 8(a) was obtained by suppressing transition (r = 0). Grid-point spacing in 
the x-direction was varied from Ax = 0.03 to 1.22 cm, and convergence to the accuracy of 
the experimental data was obtained for all values; because of the abruptness of transition, 
step sizes greater than Ax = 1.22 cm were not studied. 

Comparisons of the numerical results with typical experimental velocity profiles 
are presented in figure 8(b) for the flat-plate model and figure 8(c) for the tunnel diffuser 
floor. The agreement is seen to be very good. 



(a) Comparisons with skin- friction coefficient data. 
Figure 8.- High Reynolds number turbulent flow. 
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Tripped Turbulent Boundary Layers 


In most supersonic and hypersonic wind-tunnel facilities, it is necessary to trip the 
laminar boundary layer on small-scale models in order to simulate the full-scale condi- 
tions where most of the boundary layer would be turbulent. An example of turbulent data 
obtained by tripping the laminar boundary layer is that of Coles (ref. 87). These data 
were obtained in the Jet Propulsion Laboratory 20-inch supersonic wind tunnel. The test 
model was a sharp -leading-edge flat plate. The free-stream Mach number was varied 
from 1.966 to 4.544. Test numbers 30, 20, and 22 (see p. 33 of ref. 87) were selected 
as typical test cases. For these three cases the laminar boundary layer was tripped by a 
fence located at the leading edge of the flat plate. (See fig. 40 of ref. 87.) The skin fric- 
tion was measured at three surface locations with a floating-element balance. Boundary- 
layer profiles were measured at x = 54.56 cm. 

The test conditions for the three cases are listed as follows: 


Coles' 
test number 

M e 

p t,oo> mnAh 2 

T t,oo> K 

T w/ T t,.o 

30 

1.982 

0.09468 

302.8 

0.8295 

20 

3.701 

.1373 

311.7 

.7152 

22 

4.554 

.3894 

307.8 

.6764 


Transition was assumed to occur near the leading edge of the plate for the numerical 
calculations Xj. ^ = 0.15 cm and to be completed (T = 1.0) at Xj. j = 0.30 cm. Twenty 
equally spaced x-solution stations were used in the region 0 § x S 0.30 cm; the x-solution 
stations were then equally spaced 0.30 cm apart over the remainder of the plate. The 
total number of grid points in the 77-direction and the grid-point spacing parameter K 
were assigned values of 301 and 1.04, respectively. The computer processing time per 
test case was approximately 3 minutes. The numerical results are compared with the 
experimental velocity and Mach number profile data in figure 9 for the three test cases. 
The agreement between the numerical results and the experimental data is very good 
for all three test cases. In particular, it should be noted that the experimental skin- 
friction coefficients (see figs. 9(a), 9(c), and 9(e)) were predicted to within 1 percent, 
which is well within the accuracy range of 2 percent as quoted for the data in refer- 
ence 87. The slight disagreement between the numerical results and the Mach number 
profiles for increasing free-stream Mach number may be due in part to eddy conductivity 
effects caused by the simple Np r ^ model utilized in the calculations. 
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(b) Mach number profile for = I.982. 

Figure Comparisons with experimental data for tripped turbulent boundary layers. 



Exp. data; M e = 3.701 


Present solution 
C T0.00I62 (exp.) 

' e L 0.00 159 (predicted) 




G Exp. data; M e = 4.554 


Present solution 

c f O.OOI22 (exp.) 

,e L0.00I20 (predicted) 
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(e) Velocity profile and skin- friction coefficient for 


M 

< 



(f) Mach number profile for = 4.554. 


Figure 9* - Concluded. 



Laminar Flow With Mass Injection 

An example of data for laminar flow with transpiration cooling is presented by 
Marvin and Akin (ref. 88). These data were obtained over a range of injection rates for 
a sharp-tipped 5° cone. The cone was solid for x < 9.53 cm; the remainder of the cone 
was porous. The test conditions were as follows: 

Moo = 7-4 

P t>co = 4.137 N/m 2 

T t>co = 833.3 K 

^L=0.38 

T t 

l 9 oo 

The air-injection parameter F ranged from a minimum value of zero (no injection) to 
a maximum value of 1.3903 x 10"^. 

A comparison of the heating rate at the wall normalized by the heating rate at the 
wall just prior to the ramp injection (x = 9.53 cm) is presented in figure 10(a). The mass 
injection distribution was not uniform, as can be seen from the plot of the actual injection 
distribution normalized by the ideal injection rate. The nonuniform distribution of mass 
injection was utilized in the numerical solutions. The agreement between the numerical 
results and the experimental heat -transfer data is excellent over the entire injection 
range. Marvin and Sheaffer (ref. 89) have also obtained numerical solutions for this test 
case which are in very good agreement with the experimental data. 

This particular case is an example where variable step size must be used in the 
region of the ramp injection (x/L = 1.0). For this case the Ax values (grid-point 
spacing in the x-direction) were constant and equal to 0.3 cm up to x/L = 0.99 at which 
point the step size was decreased to a value of 0.03 cm through the ramp injection region. 
The step size was then progressively increased up to a value of 0.3 cm at x/L = 1.5. 

The flow was laminar; consequently, K was set to unity and 101 equally spaced grid 
points were used in the ^-direction. The digital computer processing time per test case 
was approximately 2 minutes. 

The effect of mass injection on the vorticity Reynolds number is shown in fig- 
ures 10(b) and 10(c). In figure 10(b) the maximum value of the vorticity Reynolds num- 
ber x max is seen to increase with increasing x. In particular, the y/6 value at 
which Xmax occurre< i is in excellent agreement with the location of the critical-layer 
position. (See ref. 72.) The effect of mass injection on the vorticity Reynolds number 
distribution is presented in figure 10(c). Increasing mass injection increases the value 
of Xmax a * a gi yen x-station as well as moves the location at which the maximum occurs 
( y / 5 )x max t° warc * the outer edge of the boundary layer. 
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(a) Comparisons with heat- transfer data. 

Figure 10.- Hypersonic laminar boundary- layer flow with mass injection. 



(b) Vorticity Reynolds number distribution for zero mass injection. 



(c) Vorticity Reynolds number distribution for variable mass injection. 

Figure 10.- Concluded. 
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Laminar Blunt Body Flow 


The numerical procedure as presented in the present paper has been applied to a 
number of blunt-body flows at supersonic and hypersonic speeds. A typical case for a 
hemispherical nose is presented in figure 11. The test conditions for this case were as 
follows: 

Moo = 10.4 

p, = 10.77 MN/m 2 
t ,00 

r n = 0.08891 m 

T+ „ = 1222 K 

LjOO 



Solution for this particular example has been obtained by Marvin and Sheaffer (ref. 89) 
and Clutter and Smith (ref. 90). The heating rate and shear stress at the wall referenced 
to the maximum heating rate and shear stress, respectively, are compared with the 
results presented in references 89 and 90 in figure 11. 

Equally spaced grid points were used in the ^-direction. Solutions were obtained 
for Ax values ranging from 0.03 cm to 0.60 cm in order to check for convergence. 



Figure 11.- Hypersonic blunt-body flow. 
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The data shown in figure 11 were calculated by using Ax = 0.15 cm although much larger 
steps could be used for engineering calculations, for example, Ax = 0.3 cm. The grid- 
point spacing in the ^-direction was held constant (K = 1) . 

Another example of hypersonic laminar flow is presented in figure 12. These 
experimental data were obtained on a spherically blunted, 25° half-angle cone. (See 
ref. 91.) The test conditions were as follows: 

M 00 = 7.95 

p. = 6.31 MN/m 2 

t,oo 

T t> „ = 783 K 



The pressure distribution used in the numerical solution was obtained by the technique 
presented in reference 92 and is compared with the experimental data in figure 12(a). 



(a) External pressure distribution. 

Figure 12.- Laminar flow over a spherically blunted 25° half-angle cone. 
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(b) Comparison with heat- transfer data. r n = 0.10} d = 1.27 cm. 

Figure 12.- Concluded. 

The numerical results for the heat-transfer coefficient distribution referenced to the 
stagnation-point value are compared with the experimental data in figure 12(b) and also 
with the local similarity method as reported in reference 91. The present results are in 
somewhat better agreement with the data than the local similarity results presented in 
reference 91. In order to assess the difference between the two predictions, the present 
method was utilized as a local similarity method by eliminating the nonsimilar terms 
from the governing equations. (See eqs. (B47) to (B49).) The results of this solution are 
also shown in figure 12(b) and indicate that the discrepancies in the present nonsimilar 
solution and the local similarity solution presented in reference 91 are not due to the 
nonsimilar terms in the governing equations. 

Nonsimilar Flow With Transverse and Longitudinal Curvature 
and Variable Pressure Gradients 

Turbulent boundary -layer data for flows where variable pressure gradients exist 
are few in number. A good example of a case where both favorable and adverse pressure 
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gradients occur as well as where transverse curvature effects are important is the data 
of Winter, Rotta, and Smith (ref. 93). The model used in the study was an axisymmetric 
configuration with continuous derivatives. (See fig. 13(a).) 


Experimental data are presented in reference 93 for the turbulent region; however, 
the solutions presented herein were obtained by starting the calculations at the tip of the 
sharp-cone forebody (x = 0). This particular configuration has received considerable 
attention over the past 2-year period. Calculations have been made by Herring and 
Mellor (ref. 94), Cebeci, Smith, and Mosinkis (ref. 83), and Beckwith (ref. 17); however, 
these solutions were all started by utilizing the experimentally measured profiles at sta- 
tion z/c = 2. The solutions presented herein are believed to be the first obtained with- 
out any dependence whatsoever on experimental profile or skin-friction data. 


The test conditions for the two cases considered are as follows: 


Moo 


Pt,oo> MN/m2 
Tt, =o, K . . 

llL 

T t ,oo 


Case 1 Case 2 

1.398 1.70 

0.0441 0.0475 

297.8 297.8 

0.976 0.971 


Z 


c 



z /c = 0^ c = 1.524 m 

r 0 /c= 0.36397 (z/c) 

z/c = O.I4I58 

r 0 /c =6.3473 (z/cf- 6.49221 (z/cf+l. 22668 (z/cf+ 0.33498 (z /c)- 0.00461 


z/c = 0.45725 

r 0 /c = - 0.236382 (z/c)t 0.60762 (z/c) 2 - 0.595521 (z/c)t0.22694 
z/C = 0.64473 

r 0 /c = -53.58258 (z/c) + 150.87806 (z/c) -I58.03866(z/cf+72.96859(z/c)-l2.4939 
z/c =0.76317 

r 0 /c = 0.236382(z/c) 3 + 0. 390789 (z/c) 2 + 0.290242 (z/c) -0.064001 
z/c = t.O 


(a) Geometry of configuration. 


Figure 13. - Comparisons with data for highly nonsimilar supersonic 
flow with transverse curvature effects. 
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The experimental Mach number distributions are presented in figure 13(b); these dis- 
tributions were used as an input to the digital computer program instead of p e . 

The initial profiles required to start the finite-difference solution were obtained by 
an exact solution of the similar boundary- layer equations (eqs. (B47) to (B49)) near x = 0. 
Transition was initiated at the x-station where x max achieved a value of 2500. The 
transition extent was then computed from equation (63). The grid-point spacing in the 
^-direction varied from a maximum Ax value of 0.3 cm to a minimum value of 0.03 cm 
in the regions of large pressure gradients. Variable grid-point spacing in the ^-direction 
was required with a K value of 1.04. Calculations were made for K values of 1.03 
and 1.05 to insure convergence. The computer processing time per test case was approx- 
imately 4 minutes. 

The numerical results are compared with the experimental data for momentum 
thickness and skin-friction coefficient distributions in figures 13(c) to 13(f). The agree- 
ment between the numerical and experimental momentum thickness and skin-friction 
coefficient distributions is good for both test cases. In particular, note the agreement 
with the minimum Cf e data point in transition (fig. 13(f)). It is also of interest to note 
that although the solutions with transverse curvature were in closest agreement with the 
0 values for x/L < 3.5, the solutions without transverse curvature agreed best with the 
C f e values in the same region. A similar trend was noted in reference 83. 


Laminar, Transitional, and Turbulent Flow on Sharp Cones 

Fischer (ref. 61) studied the effect of unit Reynolds number R e on transition on 
a sharp 10° half-angle cone. The test conditions were as follows: 


Moo = 7 

p t ^ = 1.384 to 4.199 MN/m 2 
T t ,oo * 572 K 



0.52 


The boundary -layer edge values were obtained from reference 95. The experimental 
location of transition and the extent of transition were used in the numerical calculations. 
The computer processing time per test case was approximately 2 minutes. 


Comparisons of the numerical results with the experimental Stanton number dis- 
tributions are presented in figures 14(a) to 14(f). The value of (x m ax) the 

experimental transition location is noted in each figure. The recovery factor used to 


compute the adiabatic wall temperature (see ref. 23) was as follows: 



( 75 ) 
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(b) ^ = 7j R e = 1.293 X 10 7 per meter; (>miax) cr = 2l+ 5 0 * 

Figure 14.- Comparisons with experimental Stanton number distributions 
for hypersonic flow over sharp-tipped cones. L = 2.54 cm. 




0 2 4 6 8 

x/L 


(c) Moo = 1 ; Re = 1-739 X 10 7 per meter; (x max ) c 



(d) Jk = 7; Re = 2.195 x 107 per meter; (x max ) c 
Figure 14-.- Continued. 
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(e) = 7 ; R e = 2.^54 x 107 per meter; ( X max ) cr 



Figure l4.- Continued. 



O Exp. data 
Present solution 



(h) = 8 J R e = 5-6^3 X 107 per meter? (* max ) cr 


Figure 1^-.-* Concluded. 



where r is the intermittency. (See eq. (57).) The experimental and numerical heat- 
transfer- coefficient distributions were reduced to the Stanton number distributions as 
presented in figure 14 with a common recovery factor; that is, the recovery factor 
obtained from equation (75). 


Two similar cases are presented in figures 14(g) and 14(h). These data were 
obtained at the Langley Research Center on a 10° sharp cone model. The test conditions 
for these data were as follows: 


Moo 

p, , MN/m^ 

l.OO 7 


T W 

^t,« 


Case 1 

Case 2 

8 

8 

13.93 

17.38 

809.9 

759.4 

0.40 

0.42 


The local values of M e , Rg, and (Xmax) C r are inc *i ca ted in the figures; (x ma x)cr 
ranged from 2450 to 3015. The agreement in the peak heating region is good; in par- 
ticular, note the predicted "overshoot" behavior near the end of transition and its agree- 
ment with the data. 


Laminar, Transitional, and Turbulent Flow Over a Hollow Cylinder 

O'Donnell (ref. 96) studied laminar, transitional, and turbulent boundary -layer flows 
over a hollow cylinder. Total-pressure profiles were measured at various stations along 
the cylinder. The test conditions were as follows: 

Moo = 2.41 

p t ^ = 23.7 to 406.5 kN/m2 

T t = 311.1 K 

1,00 

-^5L=0.9 

T t,00 

For this particular set of calculations, the experimental transition location was 
utilized; however, the extent of transition was calculated from equation (61). Conse- 
quently, the only variable inputs to the computer program were the specific values of the 
total pressure p^. ^ and the transition location x^. 

The velocity profile comparisons are presented in figures 15(a) to 15(e). For a 
unit Reynolds number of 2.2 X 10® per meter, the boundary layer was laminar through- 
out the measured area. The velocity profiles are similar and the agreement between the 
numerical results and experimental data is very good. (See fig. 15(a).) For unit Reynolds 
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(c) Rg = 1.89 x 10 ? per meter. (d) R e = 2.835 X lO? per meter. 

Figure 15.- Comparisons with velocity profile data and momentum, thickness for laminar 
transitional, and turbulent flow over a hollow cylinder at = 2 .^- 1 . L = 2 . 5 ^ Cl 



Momentum thickness, 9/ L 


Data Re/m 

O 2.20: 

□ 9.45 

O 18.90 

A 28.35 

A 37.73 



numbers of 9.45 x 10® per meter and 1.89 x 10 7 per meter, laminar, transitional, and 
turbulent flow occurred. (See figs. 15(b) and 15(c).) For unit Reynolds numbers of 
2.835 x 10 7 per meter and 3.773 x 10 7 per meter, the flow was turbulent as presented in 
figures 15(d) and 15(e), respectively. The experimental and calculated turbulent profiles 
are similar. Comparisons with the experimental momentum thicknesses are presented 
in figure 15(f). The agreement is good over the entire unit Reynolds number range. 

Flat-Plate Hypersonic Row 

A final example of laminar and transitional boundary-layer flow for hypersonic test 
conditions is presented in figure 16. These data were obtained by Johnson (ref. 97) on a 
sharp-leading-edge flat-plate model. The experimental test conditions were as follows: 

Moo = 7.8 

Pj = 4.226 MN/m 2 

Tf = 794.4 K 

= 0.388 

T t,oo 

The agreement between the numerical results and the experimental Stanton number dis- 
tribution is very good provided the experimental transition location is used as an input 
quantity. The transition extent was calculated from equation (63). (Note that the transi- 
tion from laminar to turbulent flow was not completed; that is, xj. f was not located on 



.3 .6 .9 1.2 

x/L 


Figure 16,- Comparison with experimental data for laminar and transitional 
flat-plate flow at i^o = 7*8. L = 30. 48 cm. 


inn i all mi | | || | 
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the model.) The total number of grid points in the ^-direction and the grid-point-spacing 
parameter K were assigned values of 201 and 1.04, respectively. The computer pro- 
cessing time was approximately 2 minutes. 

This particular test case is presented to emphasize the fact that the stability index 
(Xmax)cr is no * an invar i an f as originally suggested by Rouse (ref. 69) but is, in fact, a 
strong function of the unit Reynolds number. For all the previous test cases presented 
herein, (Xmax) C r has varied over a relatively small range, 2500 ^ (x max ) cr = 3000 
which might lead one to the erroneous conclusion that, although not invariant, (x max ) cr 

varies only slightly. However, the present flat-plate test case value of (x max ) cr = 4000 
clearly extends the range from 2500 to 4000. The numerical results obtained by assuming 
that transition would occur at the x-station where X max equaled 2500 are presented in 
figure 16. These results indicate that no assumption concerning the location of experi- 
mental transition should be made on the basis of the magnitude of the stability index as 
being an invariant or a near constant. An area of interest concerning the possibility of 
using X m ax as a transition correlation parameter would be to apply the present numer- 
ical method to particular classes of flow over a broad range of test conditions; for 
example, boundary-layer flows with nonsimilar wall boundary conditions, flows with 
pressure gradient histories, and similar flows where experimental transition data are 
available. From these systematic calculations, one could determine the possibility of 
using X m ax as a gu^e f° r predicting transition. (See ref. 23.) 

CONCLUDING REMARKS 

A system of equations which govern the mean flow properties of laminar, transi- 
tional, and turbulent compressible boundary layers for either planar or axisymmetric 
flows and a numerical method by which the system can be accurately solved has been 
presented. 

The turbulent boundary layer was treated by a two-layer concept with appropriate 
eddy viscosity models for each layer to replace the Reynolds stress terms. A specifiable 
turbulent Prandtl number replaces the turbulent heat flux term. A constant static turbu- 
lent Prandtl number equal to 0.9 was used in the present solutions; however, the numeri- 
cal procedure is completely general and any functional relation may be used. 

The mean flow properties of the transitional boundary layer were calculated by 
multiplying the eddy viscosity by an intermittency function based on the statistical produc- 
tion and growth of the turbulent spots. The intermittency function represents the fraction 
of time that any given stream wise station experiences turbulent flow or the probability at 
any given instant of time that a specific point will be engulfed in a turbulent spot. 
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The numerical method used to solve the system of equations was a three-point 
implicit finite-difference scheme for variable grid-point spacing in both spatial coordi- 
nates. The method is self starting; that is, it requires no experimental data input except 
the empirical or experimental location and extent of transition. The momentum and 
energy equations remain coupled and are solved simultaneously without iteration. The 
technique is inherently stable provided the eddy viscosity distribution is properly treated; 
that is, no constraint is placed on the grid-point spacing by a step-size stability param- 
eter such as in the case of explicit finite-difference schemes. The method is highly 
efficient with regard to flexibility, digital computer processing time, and accuracy. 

A number of test cases have been compared with experimental data for supersonic 
and hypersonic flows over planar and axisymmetric geometries. These test cases 
included laminar, transitional, and turbulent boundary-layer flows with both favorable 
and unfavorable pressure gradient histories. Mass flux at the wall and transverse curva- 
ture effects were also considered. Excellent agreement between the numerical results 
and experimental data was obtained for all test cases. However, low Reynolds number 
and high Mach number turbulent flows, supersonic flows with large favorable or unfavor- 
able pressure gradient histories, and flows with large longitudinal curvature have not 
been considered. The omission of these particular classes of flow was deliberate and 
was done for two reasons: (1) to minimize the empiricism contained in the eddy viscosity 

and eddy conductivity models and (2) the lack of accurate, well-documented data for these 
classes of flow. Advances in the state of the art of formulating turbulent transport 
models which contain low Reynolds number effects, pressure gradient effects, and lon- 
gitudinal curvature effects can be directly applied with little or no effort to the technique 
presented herein as they become available. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., June 21, 1971. 
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APPENDIX A 


DIFFERENCE RELATIONS 


Three-point implicit difference relations are used to reduce the transformed 
momentum and energy equations (eqs. (26) and (27)) to finite-difference form. It is 
assumed that all data are known at the solution stations m - 1 and m. (See fig. 7.) 

One then wishes to obtain the unknown quantities at the grid points for the m + 1 station. 
In the following development, the notations G and H are utilized to represent any 
typical variable. 


Taylor series expansions are first written about the unknown grid point (m + 1, n) 
in the ^-direction as follows: 


a t 2 


G m,n - G m+l,n - *Z 2 { G i) m+ l, n + ^“( G ^) 


m+l,n 


A£“ 


/m+l,n 


(Ala) 


and 


G m-l,n = G m+l,n - ( a «l + At -2)( G l) nHl 


( A h * 


( G «)m + l,n 


.M, c v 

6 \ ^wm+l,n 


(Alb) 


where subscript notation has been utilized to denote differentiation; that is, = (3G/9|), 
and so forth. 

Equations (Ala) and (Alb) can be solved to yield 


lBC\ _ X l G m^l,n - X 2 G m,n ♦ * A h) r . 

Wm+l,n 2A| 2 6 

and 


G m+l,n “ X 4 G m,n “ x 5 G m _! 



(A2) 


(A3) 


Terms of the order Ai^ A£ 2 or smaller are neglected. This procedure produces trun- 
cation errors of the order Ai^ A£ 2 instead of Aij 2 as in reference 29 where two-point 
difference relations are used. The Xj,X 2 , . . ., X5 coefficients appearing in equa- 
tions (A2) and (A3) are defined as follows: 
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APPENDIX A - Continued 


^A£ 1 + 2A£ 2 
Xl 2 A?l + a ^2 

(A 4 ) 

A ^1 + a ^2 

X 2 - 2 

(A 5 ) 

A£ 2 

(A6) 

X3 " 2 A^(A^ + A^ 2 ) 

A ^1 +A ^2 

(A 7 ) 

54 " Aq 

A?, 

(A8) 

X5 = A|f 


Taylor series expansions are next written about the unknown grid point (m + 1 , n) 
in the 77-direction as follows: 


A?7 n A?7 n 

G m+1 ,n+l = G m+1 ,n + A77 n( G *7) m+1>n + ~{ G w) m+ i )n + TW m+ l )n + * ' * 


(A 9 a) 


^ ^^n -1 

G m+l,n-l G m+l,n A ^n-l( Gj ?) m+ i >n + 2 ( G ^^)m+l,n 6 ( G W7) m+ i jn 


+ . 


(A 9 b) 


Equations (A 9 a) and (A 9 b) can be solved to yield 

(“ 2 ) = Y l G m+l,n+l " Y 2 G m+l,n + Y 3 G m+l,n-l 

\ d7 T 7 m +l,n 

and 

(l|)m +1 n = Y 4 G m+l,n+l " Y 5 G m+l,n " Y 6°m+l,n-l ' 


, K-l - H.) 

"t z *r 


777777 


(A 10) 


A Vl 


^777777 ■*■••• 


(All) 


The Y-pYg, . . Yg coefficients appearing in equations (A 10 ) and (All) are defined as 
follows: 


Yt = 


1 A %^n + A \-l) 


(A 12 ) 
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APPENDIX A - Continued 


Y 2 = 


y 3 = 


Y 4 = 


Ar? n At^ 


A? ?n-l( Ar ?n + A 7 ?n-l) 


Arj. 


n-1 


A77 n (A77 n + A Vn _ x ) 


(A13) 

(A14) 

(A15) 


Y 5 = 


A7? n _i ' A U 


n 


A? ?n A? ?n-1 


(A16) 


y 6 = 


An 


n 


A Vl( A?7 n + A Vl) 


(A17) 


For the case of equally spaced grid points in the £ and 77 coordinates, equa- 
tions (A4) to (A8) and (A12) to (A17) reduce to the following relations: 


= 3 
X 2 = 4 

x 3 = i7 

X 4 = 2 
X 5 = 1 


(A18a) 


and 


Y l = 


1 

At?2 


~\ 


y 2 = 2 y, 

y 3 = ^i 

y 5 = o 
y 6 = y 4 


(A18b) 




where A£ and A 77 represent the spacing between the grid points in the | and 77 
coordinates, respectively. 
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APPENDIX A - Continued 


Equations (A2), (A3), (A10), and (All) can then be written for constant grid-point 
spacing as follows: 


9G \ _ 3G 

9 £/m+l,n 


m+l,n " 4G m,n + G m-l,n , A£* 

2 — — + — — Cjrtfct + . . . 

2A| 3 «« 


(A19) 


J m+l,n = 2G m,n “ G m -l,n + A ^ G H + 


(A20) 


( 9 2 g\ _ G m+l,n+l " 2G m+l,n + G m+l,n-l n 

WL + l,n~ A, 2 '12 nrm •** 


(A21) 


9G\ _ u m+l,n+l - G n+1 ,n-l 

^Ll,n = 2 A, 


(A22) 


Equations (A19) to (A22) are recognized as the standard relations for equally spaced grid 
points. (See, for example, ref. 30.) 

Quantities of the form (g -gj-j that appear in the governing equations must be 

linearized in order to obtain a system of linear difference equations. Quantities of this 
type are linearized by using equations (A2) and (A3); that is, 

f L, n - (*«W - * °K *H) 


2 A^ 2 


(A23) 


The procedure used to linearize nonlinear products such as (9G/9?7)(9H/9?7) is 
the same as that used by Flvigge-Lotz and Blottner (ref. 29) and is as follows: 


9G\/9H 

9?7 A 07 ? 


9G 


/Jm+l,n \ dr l An ,n\ d7 l ) m +l ; n \ 97 7/m,n( 8 7?)m,n { d ^m,J\ dr l) m +l,n 


/9G\ /9H\ 


(A24) 


where the terms ( — ) and (— | are evaluated from equation (All), but at the 

\ dr l)m,n \ 9r ?/m,n 

known station m. By equating G to H in equation (A24), the linearized form for 
quantities of the type is obtained; that is, 
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APPENDIX A - Concluded 


/3G\2 = (8G\ J 

\ 9T ?/m+l,n \ 97? 'm,n \ 


(*i) 

^ / m+l,n 



(A25) 


w here (—1 is obtained from equation (A22). 

Wm+l,n 

The preceding relations for the difference quotients produce linear difference equa- 
tions (eqs. (66) and (67)) when substituted into the governing differential equations for the 
conservation of momentum (eq. (26)) and energy (eq. (27)), respectively, since terms of 
order (A|) are neglected. 
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APPENDIX B 


COEFFICIENTS FOR DIFFERENCE EQUATIONS 


Equations ( 66 ) and (67) are the difference equations used to represent the partial 
differential equations for the conservation of momentum and energy, respectively. These 
equations are repeated for convenience as follows: 


A^F i 
"n m+l,n- 


r(1)tt . p(^)w . 

B n F m+l,n + m F m+l,n+l + n °m+l,n-l 


+ E 


(■*■)© 1 + F^0 1 * = 

n m+l,n n m+l,n+l n 


(Bl) 


/v(2)p + 

A n *m+l,n-l + 


B 


( 2 ) 

n 


■<( 2) 1 


>(2), 


F m + l,n + C n‘' F m + l,„ + l +D n‘' e m + l,n. 


+ E 


(2)q . p(2)q _ q(2) 

n rn+l,n + n U m+l,n+l n 


(B2) 


These equations are obtained from equations (26) and (27) and the difference quotients 
presented in appendix A. The coefficients A^, B^, and so forth, in equations (Bl) 

and (B2) are functions of known quantities evaluated at stations m and m - 1. (See 
fig. 7.) Therefore, equations (Bl) and (B2) can be solved simultaneously without iterative 
procedures. The coefficients are as follows: 

= Y 3 H 3 - Y 6 H n (B3) 


B n } = X 1 H 1 - Y 2 H 3 - Y 5 H U + H 5 


(B4) 


= Yi h 3 + y 4 h u 

(B5) 

! n ’ = - y 6 H 4 F Y 

(B 6 ) 

i 1 * = — + H, 

n Y 6 n 6 

(B7) 

,( 1 ) _ _ Y 4 ( 1 ) 

n " y 6 n 

(B 8 ) 

= H l F m2 + H 4 t Y F Y 

(B9) 
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I APPEl'DDIX B - Continued 

f 


Ai 2) = -2Y 6 H 8 F Y 

(BIO) 

b ( 2).!5 a ( 2 ) 
B n "Y 6 A n 

(BU) 

c ( 2 )--Ii A (2 ) 

“ - y 6 " 

(B12) 

D ® - Y 3 H 10 - Y 6 H 12 

(B13) 

E n ’ = X 1 H 1 - Y 2 h 10 - Y 5 h 12 

(B14) 

F i, 2) = Y 1 H 10 + Y 4 H 12 

(B15) 

G il 2> - H l T m2 + H B< F y) 2 + H 9 (T Y ) 2 

(B16) 


The coefficients Y V Y 2 , . . Y 6 , X 1? . . X 5 , etc., are functions of the grid- 

point spacing and are defined in equations (A12) to (A17) and (A4) to (A8), respectively. 
The coefficients Hi, H 2 , . . . are defined as follows: 


■pi _ t t? FT 

H 1 ^ 

(B17) 

h 2 = V m l “ L ml( E ml c ml + E ml c ml) 

(B18) 

H 3 = " E ml L ml C ml 

(B19) 

tt _ tt L ml 

H4 ' H3 ^I 

(B20) 

h 5 = ^m+l F ml 

(B21) 

h 6 = ~^m+l 

(B22) 

h 7 = V ml ‘ L ml( E ml C ml + E ml C ml) 

(B23) 

H 8 = ' Q! m+l L ml E ml C ml 

(B24) 
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APPENDIX B - Continued 


H 9 " “ E ml L ml C ml 

(B25) 

tt _ tt L ml 
H 10 ~ h 9 t . 

(B26) 

L ml 


H n = H 2 + H 4 T y 

(B27) 

H 12 = H 7 + 2H 9 T y 

(B28) 


The undefined quantities appearing in equations (B17) to (B28) are defined as follows: 


F ml = X 4 F m,n " X 5 F m-l,n 

(B29) 

^ml = x 4 0 m,n “ ^5®m-l,n 

(B30) 

^ml = ^4^m,n ” ^5^m-l,n 

(B31) 

F m2 = X 2 F m,n “ X 3 F m-l,n 

(B32) 

T m2 = X 2 0 m,n- x 3 0 m -l, n 

(B33) 


1 + 


e/ m+l 


Lml = ^ T ml ♦ (-1 %-\ 


T e/m+l 


l' = 


T 

’(£) - T ml 

L ml 

v e Wi 

2T ml 

(f) + Tml 


Ae/ m+ i J 


■^ml (^ av ) m +l,n 
v - ( € av )m+l ,n 


Np r 


(Air only) 


(B34) 


(Air only) 


(B35) 


(See eq. (74)) 


(B36) 

(B37) 
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E Y = Y 4 i m,n+1 " Y 5 i m,n ' y 6 ? m,n-l 
E Y = Y 4 g m,n+1 ' Y 5 g m,n ' Y 6 g m,n-1 
F Y = Y 4 F m,n+l " Y 5 F m,n “ YgF m>n _4 
Ty = Y40 m>n+ 1 - Yg© m>n - Yg0 mjI1 _l 


^m+1 


= (H du e\ 

\ u e L+l 


a m+l “ 


e 'm+1 


(See eq. (All)) (B38) 

(B39) 

(B40) 

(B41) 

(See eq. (28)) (B42) 

(B43) 


The transverse curvature terms are contained in the quantities C m i and 
which appear explicitly in the H2, Hg, H^, Hg, and Hg coefficients. The transverse 
curvature term in the transformed plane (see eq. (30)) may be written as follows: 

t 21 - 1 + 2 «wW2i ± f e d„ (B44) 

P e u e do 

where t represents the ratio r/r Q and is a known quantity for the N-l grid points 
at station m - 1 and m. Then, the extrapolated values at m + 1, n are obtained as 
follows where the parameter C is used to represent t^: 

c ml = x 4 C m,n " x 5 C m-l,n ( B45 ) 


C ml “ Y 4 C m,n+l 


~ Y 5 C m,n 


" Y 6 C m,n-l 


(B46) 


Two quantities (symbols) as of now remain undefined. These are the code symbols 
FT and W which appear in equations (B17) and (B44), respectively. The code sym- 
bol W appearing in equation (B44) is used to either retain or neglect the transverse 
curvature terms for axisymmetric flows; that is, W = 1 or 0, respectively. For planar 
flows, the transverse curvature term does not appear since j equals 0. The code sym- 
bol FT (flow type) appearing in equation (B17) is used to either retain or neglect the 
nonsimilar terms in the governing differential equations; that is, FT = 1 or 0, respec- 
tively. If FT is assigned a value of unity, the solution to the nonsimilar equations 
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(eqs. (25) to (27)) is obtained. If FT is assigned a value of zero, the locally similar 
solution is obtained; that is, the following system of equations are solved: 

Continuity: 

— + F = 0 (B47) 

9 7 ] 

Momentum: 


v - 4-(t 2 he + /3(F 2 - ©) = 0 

9 T) dt] J ' 

Energy: 



JL(t 2 h -J- 8 ®\ 

*n\ N Pr 37? / 


- otlftl 



= 0 


(B48) 


(B49) 


The governing equations for the locally similar system are obtained from equa- 
tions (25) to (27) by neglecting derivatives of the dependent variables F, 0, and V 
with respect to the streamwise coordinate £. The capability of obtaining locally similar 
solutions is desirable in that for a given test case the locally similar and completely 
nonsimilar solutions can be obtained for the identical program inputs and numerical pro- 
cedures. Consequently, the effects of the nonsimilar terms on the boundary-layer char- 
acteristics for a particular case can be determined by direct comparison of the results 
obtained for solutions for FT = 1 and FT = 0, Respectively . 
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